# --- Setup ---
rm(list = ls())
setwd("/Users/John/Dropbox/")

# --- Load Required Packages ---
library(readr)
library(interflex)
library(tidyverse)
library(ggplot2)
library(cowplot)

# --- Load Data ---

df <- as.data.frame(read_csv("JOP_Replication_Materials/data/processed/final_dataset.csv"))

policies <- read_excel("JOP_Replication_Materials/appendix/data/cn_govt_policies_count.xlsx",
                       sheet = "policies") %>%
  pivot_longer(cols = 2:3, names_to = "unit", values_to = "count") %>%
  arrange(unit)

# --- Separate State Council and NDRC policies ---
sc <- policies %>%
  filter(unit == "State Council") %>%
  mutate(log_sc = log(count)) %>%
  select(year = Year, log_sc)

ndrc <- policies %>%
  filter(unit == "NDRC") %>%
  mutate(log_ndrc = log(count)) %>%
  select(year = Year, log_ndrc)

# --- Load the revenue data and process it ---
revenue <- read_excel("JOP_Replication_Materials/appendix/data/cn_govt_policies_count.xlsx",
                      sheet = "revenue") %>%
  filter(Year <= 2015) %>%
  mutate(log_rev = log(Central)) %>%
  select(year = Year, log_rev)

# --- Merge all datasets with the main dataframe ---
df <- df %>%
  left_join(sc, by = "year") %>%
  left_join(ndrc, by = "year") %>%
  left_join(revenue, by = "year")

# --- Interflex Model ---
intflex <- df %>%
  dplyr::select(isic, isic2, year, strategic, sqrtta, isic_year, median_share, 
                med_hhi_isic2, med_soe_isic2, log_ndrc, log_sc, log_rev) %>%
  mutate(Strategic = strategic)

intflex <- as.data.frame(intflex)

bin.ndrc <- interflex(estimator = "binning", Y = "isic_year", D = "Strategic", X = "log_ndrc", 
                      data = intflex, vcov.type = "robust",
                      na.rm = TRUE, main = "Marginal Effect of Strategic on Tech Absorption", Ylabel = "Tech Absorption Policies", 
                      Xlabel = "Total NDRC Policies (Logged)",
                      cex.main = 1, ncols=1, theme.bw = TRUE, bin.labs = F)

plot(bin.ndrc$figure)

bin.sc <- interflex(estimator = "binning", Y = "isic_year", D = "Strategic", X = "log_sc", 
                    data = intflex, vcov.type = "robust",
                    na.rm = TRUE, main = "Marginal Effect of Strategic on Tech Absorption", Ylabel = "Tech Absorption Policies", 
                    Xlabel = "Total State Council Policies (Logged)", method = "linear",
                    cex.main = 1, ncols=1, theme.bw = TRUE, nbins = 3, bin.labs = F)

plot(bin.sc$figure)


bin.rev <- interflex(estimator = "binning", Y = "isic_year", D = "Strategic", X = "log_rev", 
                     data = intflex, vcov.type = "robust",
                     na.rm = TRUE, main = "Marginal Effect of Strategic on Tech Absorption", Ylabel = "Tech Absorption Policies", 
                     Xlabel = "Central State Fiscal Revenue (Logged)", method = "linear",
                     cex.main = 1, ncols=1, theme.bw = TRUE, nbins = 3, bin.labs = F)

plot(bin.rev$figure)

# --- Extract the figures from the interflex models ---
fig_bin_p <- bin.ndrc$figure  # NDRC plot
fig_bin_l <- bin.sc$figure  # State Council plot

# --- Combine the plots side-by-side using plot_grid from cowplot ---
combined_figure <- plot_grid(fig_bin_p, fig_bin_l, labels = NULL, nrow = 1, align = "v")

# --- Save the combined figure to a PDF ---
plot_path <-"JOP_Replication_Materials/appendix/output/C_figure_12.pdf"
ggsave(plot_path, combined_figure, width = 14, height = 6)

browseURL(plot_path)

# Save the revenue figure to a PDF
plot_path <-"JOP_Replication_Materials/appendix/output/C_figure_13.pdf"
ggsave(plot_path, bin.rev$figure, width = 7, height = 6)

browseURL(plot_path)
